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ABSTRACT 

By  using  conjugate  directions  a  method  for  solving  convex  quadratic 
programming  problems  is  developed.  The  algorithm  generates  a  sequence  of 
feasible  solutions  and  terminates  after  a  finite  number  of  iterations. 
Extensions  of  the  algorithm  for  nonconvex  and  large  structured  quadratic 
programming  problems  are  discussed. 


AMS(MOS)  Subject  Classifications  90C20,  90C25. 
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SIGNIFICANCE  AND  EXPLANATION 


The  quadratic  programming  problem  is  the  following:  Given  n  x  1 
vectors  c,  a1,...,am,  numbers  b.,,...,bm  and  an  n  x  n  matrix  C,  find 
an  n  x  1  vector  x  which  minimizes  the  quadratic  function 

c'x  +  ^  x'Cx  , 

subject  to  the  inequality  constraints 

a! x  <  b.  ,  i  =  1 , . . . ,m. 
i  i 

If  C  is  the  n  x  n  zero  matrix,  then  the  quadratic  programming  problem 
reduces  to  the  linear  programming  problem. 

In  recent  years  quadratic  programming  has  become  an  important  tool  in 
optimization.  It  has  wide  applications  in  areas  such  as  statistics, 
structural  engineering,  economics  and  portfolio  analysis. 

The  contribution  of  this  work  is  an  algorithm  which  solves  the  quadratic 
programming  problem  in  a  finite  number  of  steps.  Furthermore,  an  extension  of 
the  algorithm  is  given  which  can  be  used  to  solve  large  structured  quadratic 
programming  problems . 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


A  QUADRATIC  PROGRAMMING  ALGORITHM 
M.  J.  Best*  and  K.  Ritter** 

1.  Introduction 

We  consider  the  quadratic  programming  problem 

min  {c. '  x  +  1  Cxi a.!x  <  ,  i  =  1, . . .  ,m}  ,  (1) 

where  c,  x,  a^,...,am  are  n-  vectors,  C  is  an  (n,n)  symmetric  matrix 

and  b,,...,b  are  scalars.  Prime  is  used  to  denote  transposition. 

1  m 

Let  F(x)  =  c'x  +  -^x'Cx  denote  the  objective  function  for  (1)  and 
R  =  {x  |  a  !x  <  b  . ,  i  =  1, . . .  ,m} 

denote  the  feasible  region,  x*  satisfies  the  Karush  -  Kuhn  -  Tucker 
conditions  for  (1)  if  there  are  numbers  u,,...,um  satisfying 

x*  e  R  , 

-c  -  Cx*  =  u-a,  +  . . .  +  u  a  ,  u  >.  0  , 

11  mm 

ui (ajx*  -  b.. )  =  0,  i  =  1 . m  . 

If  C  is  positive  semi  -  definite,  these  conditions  are  both  necessary 
and  sufficient  for  x*  to  be  a  global  minimizer  [1].  Let 

I(x*)  =  {i|ajx*  =  b^,  1  <  i  <  m}  . 

If  C  is  indefinite  and  in  addition  to  the  Karush  -  Kuhn  -  Tucker  conditions 
x*  also  satisfies 

u.j  >  0  ,  all  i  €  I(x*)  , 
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s'Cs  >  0  for  all  s  *  0,  with  a^s  =  0,  all  i  €  I(x*) 
then  x*  is  a  strong  local  minimizer  for  (1). 

We  first  consider  the  case  when  C  is  positive  semi  -definite  and 
present  an  algorithm  for  the  solution  of  (1)  which,  in  a  finite  number 
of  steps,  determines  either  an  optimal  solution  or  that  the  problem 
is  unbounded  from  below.  The  algorithm  is  based  on  a  new  updating 
procedure  for  conjugate  directions  when  the  set  of  active  constraints 
is  changed.  A  general  description  of  the  method  is  given  in  Section  2 
and  a  detailed  formulation  is  given  in  Section  3. 

In  Section  4,  we  discuss  the  case  of  C  being  indefinite.  Provided 
suitable  initial  data  is  used,  the  algorithm  of  Section  3  will 
determine  a  local  minimizer  of  (1)  in  a  finite  number  of  steps.  A 
procedure  is  then  given  which  will  ensure  that  the  initial  data 
requirement  will  be  met. 

In  Section  5,  we  consider  the  structured  quadratic  programming  problem: 
Minimize 

F1(x1)+F2(x2)  +  ...+Fp(xp)+F0(y)  (2a) 

subject  to 

avxi  +  b\^  -  *  v  =  mi_1  +  l,...,mi  (2b) 

for  i  =  1 , . . .  ,p  and 

b^  <  3v;  v  =  mp  +  1, . . .  ,m  .  (2c) 

Since  (2)  is  a  special  case  of  (1),  it  can  be  solved  using  the  algorithms 
presented  in  Sections  3  and  4.  However,  the  number  of  variables  may  be 
large  and  the  computational  expense  high.  Therefore,  it  is  appropriate 
to  develop  an  algorithm  which  takes  advantage  of  the  structure  of  (2). 
This  is  done  in  Section  5. 
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2.  General  Description  of  the  Algorithm 

A  /\  /s 

x  is  a  quasi  -stationary  point  for  (1)  if  x  f  R  and  x  is  an  optimal 
solution  for 


min{F(x)|ajx  =  b.. ,  all  i  e  I ( x )}  . 


To  ensure  finite  termination,  the  algorithm  determines  a  sequence  of 
quasi  -  stationary  points  having  decreasing  objective  function  values. 

Let  x.  denote  the  iterate  at  iteration  j.  Suppose  x.  is  not  a  quasi  - 
J  J 

stationary  point.  Assume  for  simplicity  that  the  first  q  constraints 
are  active  at  x..  We  begin  by  looking  for  a  quasi  -stationary  point 

J 

at  which  only  the  first  q  constraints  are  active.  If  such  a  point 
exists,  it  is  an  optimal  solution  for 


1  . 


minfc'x  +  ^x'CxIajx  =  i  =  1 . q)  . 


(3) 


We 


write  this  point  in  the  form  x-  -  s.,  where  s.  is  to  be  thought  of 

J  J  J 


as  a  search  direction,  s.  may  be  determined  as  follows.  Let 

J 


Dj  =  [ar-:-’aq’  Ccq+1’- •  •  ,Ccn^  ' 
and  assume  that 


(4) 


-q+1’ 


(5) 


has  columns  Cj,...,cn.  By  definition  of  the  inverse  matrix, 

cq+l’’‘’’cn  are  a  set  norTTial i zed  conjugate  directions  which  are 

orthogonal  to  the  gradients  of  the  active  constraints  (normalized  in 

the  sense  that  c^Cci  =  1,  i  =  q+1 . n).  Let  g,  =  c  +  Cx.  denote  the 

J  J 

gradient  of  F  at  x..  Define 
J 
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Because  x.  -  s.  is  optimal  for  (3),  s.  is  called  the  Newton  direction. 

*3  J  3 

Next  suppose  that  x.  is  a  quasi  -  stationary  point.  Let  D.  and  D,^  be 

J  J  J 

as  in  (4)  and  (5).  Then  gjc^  =  0,  i  =  q+l,...,n.  Let  u^ ,  i  =  l,...,q  be 

defined  by  (8).  If  u.  >  0  for  i  =  l,...,q  then  x.  satisfies  the  Karush- 

(  3 

Kuhn  -  Tucker  conditions  for  (1)  and  is  thus  an  optimal  solution.  Otherwise, 

suppose  Uq  <  0.  We  proceed  by  deleting  constraint  q  from  the  active 

set.  Suppose  first  that  c'Cc  =  0.  Then  we  set  s .  =  c  and  observe  that 

q  q  J  q 

F(Xj  -osj)  =  F<xj)  -og’.s.  . 

Since  u  =  -g'.s.  <  0,  F(x.  -  os.)  is  a  strictly  decreasing  linear 
Q  J  J  j  J 

function  of  a.  Either  x.  -as.  is  feasible  for  (1)  for  all  a  >  0  in  which 

J  J 

case  (1)  is  unbounded  from  below,  or  for  sufficiently  large  a,  some 
previously  inactive  constraint  becomes  active.  Next  suppose  that 
c'Cc  >  0.  Then  by  definition  of  the  inverse  matrix,  (c'Cc  )  *c 

H  H  H  H  H 

together  with  cq+i****«cn  form  a  set  of  normalized  conjugate  directions 
which  are  orthogonal  to  a^,...,aq  From  (6),  the  Newton  direction  is 

Sj  *  (9jCq)(cqCcqrlcq  ' 

In  fact,  it  is  more  convenient  to  use  s.  =  c  ,  which  is  parallel  to 

the  Newton  direction  and  account  for  the  scalar  (g'.c  )(c'Cc  )  1  in  the 

qM  q  q' 

stepsize  calculation.  Thus  when  either  c'Cc  =  0,  or,  c'Cc  >  0,  it  is 

q  q  q  q 

appropriate  to  set  s.  =  c  . 

j  q 


Assuming  x.  -  s.  is  feasible  for  (1),  we  wish  to  modify  D.  so  that 
3  3  J 

+  j  and  are  related  as  in  (4)  and  (5)  but  without  constraint  q, 

An  appropriate  way  to  do  this  is  to  replace  column  q  of  D'.  with 

3 
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d  =  (c’Cc  )  '  Cc  .  Let  the  new  matrix  be  denoted  by  D!  .  and  let 

'  q  q  q  j+i 


Dj+1  =  (ci’-'-’cn]  • 


The  Sherman  -  Morri son  formula  [4]  asserts 


i 

c-  =  c.  -  -x -  c  ,  for  i  =  l,...,n,  i  *  q  . 

d 1  c  " 

q 

Because  are  conjugate  directions, 


c 1  Cc  - 


:i  =  ci  "  c^cr  cq  =  ei  •  for  1  =  q+1 . n  ; 


i.e.,  the  normalized  conjugate  direction  columns  remain  unchanged  by 
the  updating.  Furthermore,  the  Sherman  -  Morrison  formula  again  asserts 


c  =  (c'Cc  f1/2c  . 

q  q  q  q 


Defi ne 


9jSj 

ijrir  •  1f  ^Csj  > 0  • 

+  °°  ,  if  s  '.Cs  .  =  0  . 

J  J 


Then  a.  is  called  the  optimal  stepsize  and  F(x.  -  as.)  is  a  strictly 
J  J  J 

decreasing  function  of  o  for  0  <  o  <  a..  So  far  we  have  made  the 

J 

assumption  that  x.  -  a.s.  €  R  provided  s'.Cs.  >  0.  We  now  suppose  that 
J  J  J  J  J 

this  is  no  longer  the  case.  Let  denote  the  largest  value  of  a  for 

which  x.  -  os.  €  R.  It  is  usual  to  call  o.  the  maximum  feasible  stepsize. 
J  J  J  K 

A 

An  explicit  formula  for  a.  is  readily  derived: 

J 
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a : x .  -  b. 


v  .  /  i  j  ~  ui 

.  nnr,{-14|— 


all  i  =  1,. . . ,m  with  a! s .  <  ol  . 

t  J  J 


Since  F(x.  -  as.)  is  a  strictly  decreasing  function  of  a  for 
J  J 

«**<*  ...  .  /v  A 

0  <  a  <  a.,  it  is  appropriate  to  set  x.  ,  =  x.  -  o.s.  with  o.  *  min{o.,o.} 
J  J  J  J 


We  continue  by  assuming  that 
a^x.  - 


=  a.  <  a..  Let  l  be  such  that 
J  J 


Then  constraint  £,  which  was  inactive  at  x.,  becomes  active  at  x. 

J  J+l 

An  obvious  way  to  proceed  is  to  obtain  D'  .  from  D'.  by  replacing  the 

J  *  J 

q-th  column  with  a^,.  If  a^  is  orthogonal  to  the  last  (n-q)  columns 
of  D.*,  then  it  follows  from  the  Sherman  -  Morrison  formula  that  these 

J 

conjugate  direction  columns  will  be  unchanged  by  the  update.  However, 
there  is  no  reason  to  expect  orthogonality  and  the  updating  would  then 
destroy  the  conjugate  directions. 

In  Lemma  1,  we  introduce  an  updating  procedure  which  circumvents  this 
difficulty.  We  motivate  it  as  follows.  Continuing  the  above  discussion, 


suppose 


a£cq+l  *  0  ' 

ak.  =  0  ,  i  =  q+2, . . .  ,n  ; 


i.e.,  a^  is  in  fact  orthogonal  to  all  but  one  of  the  conjugate 
directions.  Suppose  D!+j  is  obtained  from  by  replacing  column  q+1 
with  a^.  With 


Dj  +  1  ~  [al’ ' * • ,aq,a£’Ccq+2, . . . ,Ccnl  , 


the  Sherman  -  Morri son  formula  asserts  that 


Dj  +  1  =  [cl”",Cq,Cq+rCq+2’---’Cn]  • 


>-l 


the  point  being  that  the  last  n-q-1  columns  of  form  a  set  of 
normalized  conjugate  directions  which  are  orthogonal  to  the  gradients 
of  the  active  constraints  as  well  as  a^.  Although  constraint  q  has 
become  inactive,  its  gradient  is  still  the  q-th  column  of  Dj+j.  The 
situation  is  identical  to  that  when  x^  is  a  quasi  -stationary  point 
and  constraint  q  is  to  be  dropped.  As  previously  discussed,  we 

A 

continue  by  setting  s.  .  =  c 

J  s  j  +  1  q+1 


Of  course,  there  is  no  reason  to  expect  c^+^,...,cn  to  satisfy  (9). 
The  critical  idea  of  this  section  is  to  replace  cq+i>*--*cn  with  a 

A  A 

new  set  of  conjugate  directions  c  +j,...,cn  which  do  satisfy  (9). 

The  construction  procedure  is  based  on  the  following  lemma.  Note 
that  although  we  assume  in  this  section  that  C  is  both  symmetric  and 
positive  semi  -  definite,  the  lemma  requires  only  the  symmetry 
assumption . 


Let  e  denote  the  v-th  unit  vector, 
v 

Lemma  1 

Let  P  =  [pj,...,p^]  be  an  (n,k) -matrix  satisfying  P'CP  =  I  and  let 
d  be  an  n  -  vector  satisfying  P'd  *  0.  Let  v  be  any  integer  with 
1  <  v  <  k.  Define 
u  =  P'd  , 

r  /ufu  ,  i  f  p '  d  <  0  , 

e2  -  {  .  _ 

■ - /u Tu  ,  otherwise  , 
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0J  =  [ 2 ( u * u  -  62p;d)]‘1/2  , 

W  =  0l(02ev  "  U)  ’ 

P  =  P(I  -  2ww')  =  Iplt...,pk3  . 

Then 

A 

a)  p.  =  p.  -  i  =  where  p  =  WjPj +  ...  + wkpk 

b)  P‘ CP  =  I, 

A 

c)  d'p^.  =0,  i  =  l,...,k  and  i  *  v, 

A  A 

d)  span{p^,...,pk}  =  span{plf...,pk) . 

Proof: 

First  note  that  e^d  -  ®  so  that  9]  is  we^  -  defined.  Let  Q  =  I  - 
Then  Q  is  a  Householder  -  matrix  with  (see  e.g.  [2]) 


Q’Q 

-  QQ  =  I 

(10) 

Qu  = 

Hull  e  ^  . 

(11) 

a)  By  definition  of  P, 

P  =  P  -  2(Pw)w'  =  P  -  2pw '  , 

i  .e. 

Pi  =  Pi  -  2wip,  i  =  1, . . .  ,k. 

b)  By  (10), 

P'CP  =  Q'(P'CP)Q  =  Q'Q  =  I. 

c)  Using  (11)  we  have 

d'P  =  d’PQ  =  u'Q  =  Hu  He  , 
from  which  the  assertion  follows. 


m 


*■9 


2ww ' . 


1 


(17) 


VF . (x^ )  €  span  {a  I  v  £  J  .  .} 
i x  i '  K  v  1  j 


C.c  .  =0  for  every  v  with  a  .  =  0  (18) 

i  vi  vi 

The  first  condition  is  satisfied  if  is  a  quasi  -stationary  point  for 

(16).  The  second  condition  can  be  imposed  without  loss  of  generality. 

Indeed,  since  C.  is  a  positive  semi  -  definite  matrix  we  have  C.c  •  *  0 
1  i  vi 

if  and  only  if  c'  .C.c  .  >  0.  In  this  case  D7*  and  J.  -  can  be  updated 
J  vi  i  vi  ij  U  v 

as  in  Step  3  1  of  the  algorithm  resulting  in  an  =  -1. 


For  i  =  1 . p  define 


M.  .  =  7  c  .b‘ 
ij  L  vi  a  - 
vi 

where  the  summation  is  over  all  v  such  that  a  .•  £  1.  Then 


M! .a  =  b  for  all  v  €  J. . 
ij  v  v  lj 


and,  for  every  y, 


av(xi+Mij(yJ"y))  +  b>  =  Bv»  v€  Jij 
vFi^xi  +  Mij(yJ  -y))  e  span{avlv  £  J^.}  . 


In  order  to  verify  (22)  observe  that 


VF.(x<?  +M.  .(yJ  -y)  =  VF.(x^)  +  C.M.  .  (yJ  -y) 
i v  i  ij'-7  J '  i '  v  i  u w  J ' 


and  by  (17),  (18),  and  the  properties  of  D7 


c' .VF.(xJ)  =  0  and  c' .C.M.  .  =  0 
vi  1  '  l '  vi  1  ij 


for  all  v  such  that  a  •  <  0. 

vi 
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5.  Decomposition 


In  this  section  we  develop  a  decomposition  method  for  the  problem  (2), 
which  is  a  generalization  of  the  linear  problem  studied  by  Rosen  in  [3]. 

For  i  =  l,...,p  we  assume  that  xi  is  an  n..  -  vector,  y  is  an  n  -  vector  and 

F  i  (xi )  =  cjxj  +  -jx!  Cx^ ,  F0(y)  =  c'y  +  C  y 


are  convex  functions. 


In  order  to  avoid  some  technical  difficulties  we  assume  that  for  each 
feasible  solution  to  (2)  the  gradients  of  the  active  constraints  are 
linearly  independent. 

For  fixed  yJ ,  (2)  can  be  partitioned  into  p  subproblems  of  the  form 
min{F.(xi)|a^xi  <;  Bv  -  b^j,  v  =  m.^  +  1 . m.)  (16) 


which  can  be  solved  by  the  algorithm  described  in  Section  3.  If  the 
feasible  set  of  (16)  has  extreme  points  and  if  F. (x. )  is  linear  we  can 
assume  that  the  optimal  solution  x^  is  an  extreme  point.  Then  the  active 
constraints  can  be  used  to  eliminate  the  x.  variables.  If  F_.(x^)  1S  a 
quadratic  function  then  an  optimal  solution  x?  to  (16)  is  in  general 
not  an  extreme  point.  In  order  to  eliminate  all  x^  variables  we  use 
appropriate  columns  of  the  matrix  DTj  associated  with  xj. 

More  generally  let  x'j,  D~.1.  =  (c^  , . . .  ,cn  ■ ) ,  =  (a^ . ap  ^  be  a 

feasible  solution  for  (16),  the  associated  matrix  and  index  set,  respectively. 


such  that 


Following  3),  it  may  be  necessary  to  return  to  1)  and  2).  This 

process  may  be  repeated  several  times  if  necessary.  It  must  terminate 

after  a  finite  number  of  steps,  however.  At  each  application  of  1),  2) 

or  3),  the  number  of  a..'s  having  value  zero  is  decreased  by  at  least 

^  J 

one.  Furthermore,  application  of  the  algorithm  cannot  increase  the 
number  of  these  a..'s.  Therefore,  after  a  finite  number  of  steps,  a 

'  J 

Karush  -  Kuhn  -  Tucker  point  x.  must  be  obtained  with  either  a .  .  *  0 

J  '  J 

for  i  =  1, . . .  ,n,  or 

c!.Cc..  =  0  ,  for  all  i  with  a. .  =  0  ,  (14) 

*  J  I  J  1 J 

and 

CkjCCpj  =  0  '  for  each  Pair  with  =  apj  =  0  •  (15) 


Assume  x.  satisfies  the  strict  complementary  slackness  condition.  We 

J 

claim  that  if  a. .  *  0  for  i  *  l,...,n,  then  x.  is  a  strong  local 
■  j  J 

minimizer  and  if  a.  •  =  0  for  at  least  one  i  then  x.  is  a  weak  local 

'  J  J 

minimizer.  In  the  former  case,  the  argument  is  identical  to  the  proof 
of  Theorem  5.  In  the  latter  case,  let  s  be  any  n  -  vector  with  a ! s  =  0, 
all  i  e  I(x-).  From  Lenina  2c),  there  are  numbers  w.  such  that 

J  I 


s  =  y  w.c. .  + 

a. .  =  -1  1  ^ 
ij 


w.c . . 
,  1 


From  Lemma  2a),  2b),  (14)  and  (15) 


s'Cs  =  l  wf  ^  0  , 
“ij--1 1 


from  which  the  assertion  follows. 
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1)  If  akj  =  0  and  ckj.Cckj  >  0,  update  using  Step  3.1. 

2)  If  akj.  =  0  and  c^Cc^  <  0,  set  s.  =  ck^  and  proceed  with 
the  algorithm  until  a  new  Karush  -  Kuhn  -  Tucker  poirt  is 
obtained . 


Repeat  1)  until  c1.  .Cc.  .  <  0  for  all  i  with  a.  .  =  0  then  perform  2). 

^  J  ^  J  J 

Repeat  this  process  until 

c!  .Cc.  .  =  0  for  all  i  with  a. .  =  0  . 
ij  U  ij 

Additional  calculations  may  be  required  to  determine  whether  or  not  x. 

3 

is  a  local  minimizer. 

3)  If  there  are  p  and  k  with  a  .  =  a.  .  =0,  c'  .Cc  .  =  c,‘  .Cc,  .  =  0 
'  PJ  kj  pj  pj  kj  kj 

if  c'  -Cc,  .  <  0 
03  kj 

otherwise  , 

and  proceed  with  the  algorithm  using  s.,  setting 

3 

ak  =  ap  j+j  =  0  in  Step  3.2,  until  a  new  Karush  -  Kuhn  -  Tucker 
point  is  determined. 

If  s.  is  constructed  as  in  3),  then  s^Cs.  =  |2|c'.Cc.  .  <  0,  and 
J  J  3  PJ  Kj 

F(xj  '  osj>  -  F(xj>  *  °2siCsj  • 

Thus  F(x,  -  as.)  is  a  strictly  decreasing  function  of  a  for  all  o  >  0 

J  J 

and  Xj  is  not  a  local  minimizer. 


and  c'  .Cc.  .  *  0,  set 
PJ  kj 


=  /  CPJ  +  Ckj  ’ 
^  CPJ  ‘  Ckj  » 


a ! s  =  0,  all  i  €  I(x.)*  Then  from  Lemma  2c),  there  are  numbers  w.. 
^  J 

such  that 


From  Lemma  2a), 

s'Cs  =  l  (w  )2  . 

a. .  =  - 1 
1J 

At  least  one  of  the  v^'s  must  be  non -zero.  Therefore  s'Cs  >  0  and  x^ 
is  indeed  a  strong  local  minimizer.  We  have  proved 

Theorem  5 

Let  C  be  indefinite  and  let  xQ  be  an  extreme  point.  Then  the  algorithm 

terminates  in  a  finite  number  of  steps  with  either  the  information  that 

(1)  is  unbounded  from  below,  or  a  Karush  -  Kuhn  -  Tucker  point  x  ..  In  the 

J 

latter  case,  if  x.  satisfies  the  strict  complementary  slackness  condition 
J 

then  x.  is  also  a  strong  local  minimizer  for  (1). 

J 

The  assumption  that  xQ  is  an  extreme  point  is  quite  strong.  Indeed,  R 
may  not  possess  an  extreme  point.  We  now  remove  the  assumption.  Let  xQ 
be  an  arbitrary  feasible  point  and  suppose  the  algorithm  has  been  employed 
to  obtain  a  Karush  -  Kuhn  -  Tucker  point  x..  Let  DT1  and  J.  be  the  associated 

J  J  J 

data  determined  by  the  algorithm.  The  property  stated  in  Lemma  4  may  not  be 
satisfied.  Further  calculations  must  be  made  in  order  to  either  determine 
that  x-  is  a  local  minimizer  or  find  a  local  minimizer  with  a  better 

J 

objective  function  value.  Consider  the  following  two  steps. 
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is  a  quasi  -  stationary  point  and  the  lemma  again  holds  at  the  next 

quasi  -  stationary  point.  If  o.  =  o.  and  the  updating  proceeds  via 

J  J 

Step  3.2  then  =  0  and  gj+jCk  j+1  *  0.  Furthermore,  a.  *  0 

for  all  i  =  l,...,n  with  i  *  k.  At  iteration  j+1,  the  algorithm  will 
choose  in  Step  1.1,  parallel  to  column  k  of  DT*j.  Successive 
iterations  will  choose  the  search  direction,  according  to  Step  1.1, 
parallel  to  the  k-th  column  of  the  current  inverse  matrix.  Each  time 
the  updating  uses.  Step  3.2  or  3.3,  some  new  constraint  becomes  active. 

In  at  most  n  steps  therefore,  either  an  extreme  point  is  located  or 
the  optimal  stepsize  is  used.  In  the  latter  case,  the  k-th  column 
of  the  new  inverse  matrix  is  a  conjugate  direction.  In  either  case, 
the  next  iteration  is  a  quasi  -  stationary  point  for  which  the  lemma 
holds.  The  assertion  of  the  lemma  for  all  quasi  -stationary  points  now 
follows  by  induction. 

The  finite  termination  argument  of  the  previous  section  did  not  require 

that  C  be  positive  semi  -  definite.  Therefore,  with  the  non  -  degeneracy 

assumption  of  Section  3,  the  algorithm  will  terminate  in  a  finite  number 

of  steps  when  C  is  indefinite.  Now  consider  the  case  when  termination 

occurs  at  iteration  j.  Then  x.  satisfies  the  Karush  -  Kuhn  -  Tucker 

J 

conditions  with  multipliers 

u  =  -g'.c. .  ,  for  all  i  with  1  £  a. .  <  m  , 
aJ  iJ  ij 

and 


u.j  =  0  ,  otherwise  . 

Assume  that  x.  satisfies  the  strict  complementary  slackness  condition; 
J 

i.e.  u.j  >0  for  all  i  €  I(Xj).  Let  s  be  any  non  -  zero  n  -  vector  with 


4.  The  Non  -  Convex  Case 


We  now  consider  (1)  when  C  may  be  indefinite.  In  this  case,  (1)  may 
possess  several  local  minimizers  and  we  consider  the  problem  of  modifying 
the  algorithm  of  Section  3  to  obtain  one.  Because  C  is  no  longer  positive 
semi  -  definite,  it  may  occur  in  Step  2  that  sjCs^  <  0.  But  then 
F(x.  -  os.)  is  a  strictly  decreasing  function  of  o  for  all  a  k  0  and 

J  w 

setting  o.  =  +®  is  appropriate.  Suppose  we  begin  the  algorithm  with  an 
J 

extreme  point  xQ  for  (1).  We  will  show  that  with  no  further  modifications, 
the  algorithm  will  terminate  with  either  a  local  minimizer  or  the 
information  that  (1)  is  unbounded  from  below. 

Lemma  4 

Let  x^  be  an  extreme  point  and  let  x,-,  D~*  and  J.  be  determined  at  the 
o  K  J  J  J 

j-th  iteration  of  the  algorithm.  Then  for  each  quasi  -stationary  point 
x.  either  1  £  a. .  <,  m  or  a. .  *  -1  for  i  *  l,...,n. 

J  ■  J  I J 

Proof: 

Since  xq  is  an  extreme  point,  the  lemma  is  verified  for  the  first  quasi  - 

stationary  point.  We  proceed  by  induction.  Assume  that  x.  is  a  quasi  - 

J 

stationary  point  and  that  the  assertion  of  the  lemma  holds  for  it  and 
all  previous  such  points.  The  algorithm  proceeds  in  Step  1.2  by  examining 
the  Lagrange  multipliers  u  =  -g'.c..  for  each  active  constraint  a... 

a.j j  ij  ij 

Let  k  be  as  in  Step  1.2  with  glc.  .  >  0.  Then  s.  =  c.  ..  If  a.  =  a.,  then 

J  KJ  J  KJ  J  J 

Xj+1  is  a  quasi  -  stationary  point,  the  updating  proceeds  via  Step  3.1, 

cl^  j  +  j  =  -1  and  the  lemma  holds  for  the  next  quasi  -  stationary  point. 

A 

If  o,  =  o,  and  the  updating  proceeds  via  Step  3.3,  then  a.  .  .  >  1,  x.  . 

J  J  K ,  J"*"'*’  J-’* 
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for  at  most  n  iterations.  Suppose  then,  that  s.  is  constructed  from 

J 

Step  1.2.  From  Step  1.1  either  there  is  no  i  with  =  0,  or. 


q'.c.  .  =  0  ,  for  all  i  with  a.  .  =  0  . 
yJ  iJ  iJ 


(13) 


Let  D'.  =  [d..,...,d  .].  By  definition  of  the  inverse  matrix, 
J  —  J 


9j  i^9JCij^dij 


With  (12)  and  (13)  this  implies 


9-j  =  l  (g'c.)a  , 

J  1  <a-  •  < m  J  ai  j 

•  J 


which  implies  that  x.  is  a  quasi  -stationary  point. 

J 


Suppose  again  that  s.  is  constructed  from  Step  1.2.  If  a.  >  0,  then 

J  w 

F(Xj+1)  <  F ( x j )  and  the  next  quasi  -  stationary  point  determined  will 
have  an  objective  function  value  strictly  less  than  F(x.).  Then  the 

J 

associated  set  of  active  constraints  can  never  be  repeated.  Since 
there  are  finitely  many  subsets  of  the  integers  l,2,...,m,  termination 
in  a  finite  number  of  iterations  is  assured.  If  a.  =  0,  some  constraint 

J 

which  is  active,  but  not  in  the  active  set,  is  then  added  to  the 
active  set.  This  could  happen  on  several  consecutive  iterations.  With 
the  non  -  degeneracy  assumption,  however,  in  no  more  than  n  steps  all 
active  constraints  must  be  in  the  active  set.  Then  a  strictly  positive 
stepsize  must  be  obtained  and  the  previous  argument  applies. 


for  all  i  with  i  *  k  . 


ai,j+l  =  aij  ’ 

Replace  j  with  j+1  and  go  to  Step  1.1. 

The  critical  properties  of  the  matrix  DT*  are  summarized  in  Lemma  2. 

J 

They  are  easily  proved  using  Lemma  1  and  the  Sherman  -  Morrison  formula. 


Lemma  2 

Let  Dj  =  [Cy . c^.J  and 

algorithm.  Then 


a) 

cljCcij 

=  1  , 

cfjCclcj 

b) 

cijCckj 

»  0  , 

for  all 

c) 

a 1  c.  . 
aij kj 

»  0  , 

for  all 

d) 

aa,..cij 

=  1  , 

for  all 

,a  .}  be  determined  by  the 

l  lj 


=  0  for  all  i,k  with  i  *  k  and  a^.  =  = 

i,k  with  0  <  a  —  <.  m  and  =  -1  , 
i ,k  with  1  <  a. ■  ^  m  and  k  *  i  , 

*  J 

i  with  1  S  a. •  S  m  . 

*  J 


Theorem  3 

The  algorithm  terminates  in  a  finite  number  of  steps  with  either  an 
optimal  solution  for  (1)  or  the  information  that  (1)  is  unbounded  from 
below. 


Proof: 

For  each  iteration  j,  it  follows  from  Lemma  2a)  and  b)  and  Step  3  that 
g'c. .  =  0  for  all  i  with  a..  =  -1  .  (12) 

J  1  J  1  j 

If  Sj  is  constructed  by  Step  1.1  then  either  contains  an  additional 

conjugate  direction  column  (a.  =  o. )  or  the  number  of  active  constraints 

J  J 

A 

increases  by  one  (o.  =  o.).  Therefore,  Step  1.1  can  be  used  consecutively 


Step  3. 


(yj> 

i-{ 

a£ci  j  ’ 

0 

if  a. .  =  -1  , 
ij 

otherwise  . 

Set 

ui  = 

{" 

.  if  (y.)  >  0 

J 

l  +1 

,  otherwise  , 

wj  - 

II  Uj 

l|yjl|ev- 

yjirVjllyjK  -  J'j)  > 

pj = 

I 

aij 

u  • 

c.  .  .  =  c . .  -  2(w.).p.,  for  all  i  with  a..  =  -1  and  i  *  v 
1  +  1  J  *  J  '  J 

c  =  (aUc  .  -  2(w.)  p.))_1(c  .  -  2(w.)  p.)  , 
v.j+l  v  V  vj  v  jVj"  '  vj  '  jVj' 

ci,j+i  *  cij  -  <aicij)cv,j+i’  for  811  *  w1th  0  £  “ij  £  ■  • 

a.  .  ,  =  a. .  ,  for  all  i  with  i  *  v  and  i  *  k 
1  *J  +  1 


%.j+l  ■  *  • 
=  0  • 


Replace  j  with  j+1  and  go  to  Step  1.1, 


Set 

Ck,j+1 

Ci » j+1 


-1 


=  (aickj>  ckj  • 


Hi  -  [ 


CkJ  • 


for  all  i  with  i  *  k 
and  0  <  a. .  <  m  , 

*  J 


c .  . 


c, . 


for  all  i  with  a.  .  =  -1 


If  o.  «  o.  *  +«> f  print  the  message  “objective  function  is 
J  J 

^  A 

unbounded  from  below"  and  stop.  Otherwise,  set  cr.  =  min{o.  ,o.) 

J  J  J 

and  go  to  Step  3. 


Step  3:  Computation  of  x^,  and 

Set  x.+1  -  x.  -  a.s.,  gi+1  -  c  +  Cxj+r  Compute 

1 


Jj+1 


{a. 


lj+l1 
s.  If 

go  to  Step  3.2. 


,an,j+l}  and  Dj  =  [cl,j+l . c 


n,j+l 


] 


as  follows.  If  a.  =  o.,  then  go  to  Step  3.1. and  otherwise 
J  J 


Step  3.1: 


Set 
c 


k.j+1 


Ci,j+1 


Ci,j+1 


(ckjCckj)  1/2ckj 

.  rCkjCci j 

ij 


ckj. ,  for  all  i  with  i  *  k, 
and  0  £  a. .  £  m 

'  J 


c  - .  ,  for  all  i  with  a..  =  -1, 

I  J  1 J 


cxk,  j+1  =  _1  ’ 

a.  .  ,  =  a. .  ,  for  all  i  with  i  *  k  . 

*  t  J  +  1  I  J 

Replace  j  with  j+1  and  go  to  Step  1.1. 


Step  3.2:  If  a. .  2  0  for  i  =  1, . . . ,n,  then  go  to  Step  3.3. 

*  J 

Otherwise,  let  v  be  such  that 

ia£cvji =  max{ia£ciji  i a11 1  with  ajj =  -u  • 

If  ale  .  =0,  go  to  Step  3.3.  Otherwise,  for  i  =  l,...,n, 

■C  VJ 


-12- 


Step  1:  Computation  of  Search  Direction  s. 

J 

Ut  D.1  =  and  =  ta^. , . . .  ,an  .) 


Step  1.1:  If  there  is  no  i  with  ou ^  =  0,  go  to  Step  1.2. 
Otlerwise,  let  k  be  such  that 

1 9 j ck j  I  =  maxtl9jcijl  I  a11  1  with  aij  =  °J  • 

If  gjc^  =  °*  go  to  Step  1.2  and  otherwise  set 
cn»  if  9-Cn  >  0, 


5  =  J  vkj’  ”  yj"kj 


c,  • ,  otherwise, 


and  go  to  Step  2. 


Step  1.2:  If  there  is  no  i  with  1  <  a. .  £  m,  then  stop  with 

solution  x..  Otherwise,  let  k  be  such  that 
J 

9jckj  =  maxtgjc^.  |  all  i  with  1  <  <  m)  . 

!f  gjckj.  <  0,  then  stop  with  solution  Xy  Otherwise 
set  s.  =  c,  .  and  go  to  Step  2. 

J  KJ 


Step  2:  Computation  of  Stepsize  o. 

J 

Compute  siCs..  If  s'.Cs.  <  0,  then  set  o.  =  +«>. 

J  J  J  J  J 

Otherwise,  set 
5  =  9JSj 

J  sicV 


If  a^.Sj  >  0  for  i 


1, . . . ,m,  set  a.  =  +  < 
J 


Otherwise,  compute  l  and  such  that 


~  d£Xj 
O.  =  ’ 


a£s. 


r  a ! x  •  -  b. 

min  {  — ^ - 1 

1  ai  sj 


i  with  a ! s .  <  0  1 
1  J  J 


bAU  JU  mtm 


We  now  give  a  detailed  statement  of  the  algorithm.  The  initial  data 
required  is  a  feasible  point  xQ*  an  ordered  index  set  JQ  =  faio*"",aho^ 
and  an  (n,n)  matrix  D”1  =  [clQ, . . .  ,cn()] .  Letting  =  [dlQ, . . . ,dnQ] ,  the 
initial  data  must  satisfy  0  <  a-  <  m  for  i  =  l,...,n  and  for  each  i 

TO 

with  1  <  a.  <  m,  we  require  a'  x  =  b  and  d.  =  a 
10  a.  o  a-,,  10  Cl- 

10  10  10 

At  a  general  iteration  j,  the  algorithm  has  available  x.,  the  ordered 

J 

index  set  J,  =  {a,.,...,a  .}  and  the  (n,n)  matrix  DT1  =  [c. c  .]. 

J  "  J  J  llj 

Letting  =  [d. . , . . .  ,d  .] ,  for  each  i  with  1  <  a . .  <:  m,  constraint  a.. 

is  active  at  x-  and  its  gradient  is  the  i-th  column  of  D'..  All  c.  -  for 
J  J  i  J 

which  =  -1  form  a  set  of  normalized  conjugate  directions  which  are 
orthogonal  to  the  gradients  of  the  active  constraints. 

We  assume  that  each  x  €  R  is  non  -  degenerate;  i.e.,  the  gradients  of 
those  constraints  active  of  x  are  linearly  independent. 


With  initial  data  xQ,  DQ  and  Jo>  we  set  j  =  0  and  the  steps  of  the 
algorithm  are  as  follows. 
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cl)  The  result  follows  from  P  =  PQ  and  PQ  =  P(QQ)  =  P. 


Continuing  the  previous  discussion,  we  can  use 


P  =  [cq+1.---*cn]  .  d  =  a l 
and  then  apply  the  lemma  to  get 


the 


P  =  fcq+l . Cn3  • 

Part  a)  gives  c^  in  terms  of  the  c^  and  a^,  part  b)  shows  that 
c.j  are  normalized  conjugate  directions,  and  part  c)  shows  that  (9) 
is  satisfied.  Parts  b)  and  d)  show  that  if  01  is  replaced  with 

J 


then 


Dj  -  [ a ^ , . . . , a^ , CCq_^ , . . . , Cc^]  , 


It  is  possible  that  q  =  n;  i.e.,  exactly  n  constraint  are  active  at 

Xj,  or,  that  a^  is  orthogonal  to  cq+i»***»cn*  In  this  case,  we  obtain 

D'.  .  from  D'.  by  replacing  a  with  a„.  D73.  is  obtained  from  D.3  and 

J+l  J  s  q  l  J+l  J 

using  the  Sherman  -  Morrison  formula  and  the  conjugate  direction 

columns  of  DT*  remain  unchanged. 

J 

We  allow  the  algorithm  to  begin  with  an  arbitrary  feasible  point  xq. 
Suppose  constraints  1,2,..., q  are  active  at  xq.  Let 


°o  "  [ar--*’aq’dq+1’---’dn]  ’ 


where  d  d  are  any  n  -  vectors  such  that  D  is  non  -  singular. 

q+1  *  n  J  o  3 

If  q  <  n,  we  cannot  in  general  assume  that  the  last  n-q  columns  of 


Substituting 


x.  =  X;  +  M.  .  (yJ  -y) 

i  i  ij'J  J ' 


(23) 


into  F^(x^)  and  the  remaining  constraints  of  (16)  we  obtain 
Ojj(y)  :■  F,(x^  +  Mij(yj  - y)) 


and 


a*  (x'l  +M.  .(y^  -y))  <  3 
v  i  i/r-7  J ' '  \ 


or 


(b ‘  -  a ‘M.  . )y 
v  v 


This  leads  to  the  following  master  problem: 
Minimize 


Qn-(y)  l  Q;Ay)  +  F  (y) 


tru 


i=l 


subject  to 


b>  -  6v  ’  v  =  m  +l,...,m 


and 


(b '  -  a 'M .  .  )y  <  B  -  a ' (x?  +  M . .yJ ) 
v  v  v  vv  i  iy  ' 

for  i  =  1,. . .  ,p  and  all  v  e  ^  +  *»•  -  -  »m^}  such  that  v  C  J..  .. 

Let  yj  ,  DJ+1  =  (ci  J+1>...,cnJ+1),  and  Jj+J  =  {ai>J+i*'  •••an,j+l} 

be  an  optimal  solution  to  the  master  problem  and  the  associated  matrix 

and  index  set,  respectively.  Define  the  sets  I.,  I . . ,  ....  I  .  such  that 

J  ^  J  P  J 


if  r  €  Jj  +  j  and  r  £  {m..^  +  1, ... ,m..} . 

Since  gradients  of  active  constraints  of  (2)  are  assumed  to  be  linearly 
independent,  it  follows  that  the  gradients  of  the  active  constraints 
of  the  master  problem  are  linearly  independent.  Thus  we  may  assume  that 
the  gradients  of  all  constraints  which  are  active  at  y^+1  are  among 
the  columns  of  Dj+1* 

Use  (23)  to  define 

xi  1  =  xi +  MijCy3  -yj  )  »  1  =  1 . p* 

Then  it  follows  from  (21)  that 

(xf.xf1 . xj+V+1)  (24) 

is  a  quasi  -  stationary  point  for  the  problem  (2)  if  there  are  numbers 

X  ,  to.  ,  and  X.  such  that 
v  iv  l  r 

-VFi(xf1)  =  l  uivav+  l  Xirar,  i  =  l,...,p  (25) 

v€Jij  r€lij 

-VF  (yj+1)  =  £  Xb  +  f  f  7  w.  b  +  l  X.  b  1.(26) 

0  v€  I  •  v  v  i  =  lUe  J. .  lv  v  r  £  I  •  •  ir  rJ 


We  will  first  show  that  (24)  is  a  quasi  -  stationary  point  if  the  following 
condition  is  satisfied. 


ar  €  span{avlv  £  J.^}  for  all  r  £  !„,  i  =  1 . p.  (27) 


Indeed,  it  follows  from  (22)  and  (27)  that  there  are  numbers  x.  and  p 
such  that 


(32) 


(33) 


-27- 


By  (29)  and  (20)  we  have  for  i  =  l,...,p. 


-M!  .a  *  l  p  b  for  all  v  €  I. 


v£Jij 


ij 


Hence,  (26)  follows  from  (30)  -  (33). 


In  order  to  compute  a  quasi  -stationary  point  for  problem  (2)  we, 
therefore,  solve  the  subproblem  (16)  and  then  formulate  and  solve  the 
master  problem.  Finally,  we  use  Step  3.2  of  the  algorithm  to  update 
DTj  by  incorporating  as  many  gradients  a  ,  r  €  I...  into  D! .  as 
possible.  Denote  the  new  matrix  and  index  set  by  DT1^  and 
respectively.  If, 


Ji,j+1  =  Jij  f0r  1  =  1 . .  * 

then  (27)  holds  and  (24)  is  a  quasi  -stationary  point.  If 
for  at  least  one  i,  we  use  x^+1,  DT1.  .,  and  J.  ,  .  to  formulate  and 
solve  a  new  master  problem.  Since  this  iteration  strictly  increases  the 
total  number  of  positive  elements  in  the  index  sets,  a  finite  number  of 
iterations  suffices  to  generate  a  quasi  -stationary  point. 


Now  let  us  assume  that  (24)  is  a  quasi  -  stationary  point,  i.e.,  equalities 


(25)  and  (26)  hold.  We  have  seen  that  £  0  and  \^r  £  0.  Thus  (24)  is 
an  optimal  solution  if  £ 
deduce  from  (28)  -  (30)  that 


an  optimal  solution  if  £  0.  Using  the  properties  of  and  DJ+i  * 


J+l 


Tiv  =  'civVFi^xi  v€jij’  1  = 


(34) 


p  =  -c!  a  , 
rv  iv  r 


V  £  Jij’  i  =  1****’p  ’ 


J+l, 


1  ”  l  }  I  *  •  |P  • 


Thus  the  multipliers  can  be  computed  from  (31). 
For  i  =  l,...,p  determine  k.  with 


Mik.  =  min{aWv  €  3..} 

~  /V 

and  define  the  set  I.  c  {l,...,p}  such  that  i  €  I.  if  and  only  if 

J  3 

a)..  <  0.  Then  (24)  is  an  optimal  solution  to  (2)  if  and  only  if  I .  =  0. 

Ki  J 

Let  i  €  I. .  Then  further  progress  can  be  made  if  the  constraint  with  index 

vJ 

k^  is  dropped  from  the  set  of  active  constraints.  Since  (21)  shows  that 
every  constraint  with  v  €  Is  treated  like  an  equality  constraint  we 

*  *J 

have  to  delete  the  gradient  a.  from  the  matrix  D!..  If  I..  =0,  it 

k,  1J  U 

follows  from  (31)  that  x. .  =  to.,  <  0,  which  by  (34)  implies  that  x^+* 

i  1Ki  1 

is  not  an  optimal  solution  to  the  subproblem 


mintF^x^la^x.  <  Bv-b^yJ  +  1,  v  -  m.^  +  1, . . . ,m.}  . 

Thus  we  can  use  the  algorithm  of  Section  3  with  initial  data  xJ!+1, 

-1  J 
D..,  and  J..  to  compute  an  optimal  solution. 

1  J  I  J 

If  L  -  *  0,  choose  any  l  e  I.  .  and  update  D7*  as  in  Step  3.3  of  the 
J  '  J  1  j 

algorithm  by  replacing  a^  in  DI^.  with  a^. 

-1  * 

After  modifying  the  matrix  D. .  in  this  way  for  every  i  €  I .  we  define 

'  J  j 

and  solve  a  new  master  problem.  Continuing  as  described  above  we  will 
obtained  a  new  quasi  -stationary  point  which  gives  a  smaller  value  of 
the  objective  function  (2a)  than  the  previous  point  (24).  Thus  the  method 
will  terminate  after  a  finite  number  of  iterations. 
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